library(tidyHeatmap)
source('00_util_scripts/mod_bulk.R')
source('00_util_scripts/mod_bplot.R')

kegg_mva <-
  c('Hmgcr','Hmgcs1','Hmgcs2','Fdps','Mvd','Idi1','Idi2', 'Pmvk','Acat1','Acat2','Mvk') |>
  str_to_upper()

# examine interested cytokines ------
key_cytokine <- c('Il6','Ccl20','Tslp','Flt3lg','Csf2','Tnf') |>
  str_to_upper()

# xiong 19 HaCaT ---------
xio19 <- download_ncbi_counts('GSE138800')

tdb.x19 <- xio19 |>
  select(-2) |>
  set_names(c('gene','ctrl1','ctrl2','ctrl3','lUVB1','lUVB2','lUVB3','hUVB1','hUVB2','hUVB3')) |>
  pivot_longer(-1) |>
  mutate(group = str_remove(name, '.$')) |>
  tidybulk(.sample = name, .transcript = gene, .abundance = value)

tdb.x19 <- tdb.x19 |>
  preproc_bulk(group)

tdb.x19 |>
  plot_qc_bulk(value_scaled, group)

res.x19 <- tdb.x19 |>
  test_differential_abundance(~ 0 + group, contrasts = 'grouphUVB - groupctrl',
                              omit_contrast_in_colnames = TRUE)

res.x19 |>
  pivot_transcript() |>
  filter(gene %in% key_cytokine, FDR < .05)

res.x19.cyt <- res.x19 |>
  mutate(group = fct_relevel(group, 'ctrl', 'lUVB')) |>
  filter(gene %in% key_cytokine)

cusp.x19.cyt <- res.x19 |>
  pivot_transcript() |>
  filter(gene %in% key_cytokine, FDR < .05) |>
  select(gene, FDR) |>
  right_join(res.x19.cyt) |>
  customize_pvalue(y = value_scaled, facets = gene,
                   p.value = FDR, group.1 = 'ctrl', group.2 = 'hUVB')

res.x19 |>
  mutate(group = fct_relevel(group, 'ctrl', 'lUVB')) |>
  filter(gene %in% key_cytokine) |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_summary(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~gene, scales = 'free_y') +
  theme_pubr() +
  scale_color_manual(values = c('blue','orange','red'), label = c('ctrl','low-UVB','high-UVB')) +
  add_manual_pvalue(cusp.x19.cyt) +
  scale_y_continuous(expand = expansion(mult = c(.1,.1))) +
  labs(y = 'Normalized expression', title = 'RNA expression of HaCaT cells after UVB treatment',
       subtitle = 'GSE138800')

x19.mva.type <- res.x19 |>
  pivot_transcript() |>
  filter(gene %in% kegg_mva, FDR < .05) |>
  mutate(type = ifelse(logFC > 0, 'Upregulated', 'Downregulated')) |>
  select(gene, type)

res.x19 |>
  mutate(group = fct_relevel(group, 'ctrl', 'lUVB')) |>
  filter(gene %in% kegg_mva, FDR < .05) |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_summary(geom = 'errorbar', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~gene, scales = 'free_y') +
  theme_pubr() +
  scale_color_manual(values = c('blue','orange','red'))

res.x19 |>
  mutate(group = fct_relevel(group, 'ctrl', 'lUVB')) |>
  filter(gene %in% kegg_mva) |>
  as_tibble() |>
  select(gene, name, value_scaled, group) |>
  left_join(x19.mva.type) |>
  mutate(type = ifelse(is.na(type), 'NS', type)) |>
  rename(sample = name, z_score = value_scaled) |>
  group_by(group, type) |>
  heatmap(gene, sample, z_score, scale = 'row',
          column_title = 'RNA expression of HaCaT cells after UVB treatment',
          palette_value = c('blue','white','red'),
          palette_grouping = list(c('green3','grey','salmon'),
                                  c('blue','orange','red')))

# human skin bulk-seq GSE148535 ------
skin.bk <- download_ncbi_counts('GSE148535')

skin.bk

skin.bk.meta <-
  download_ncbi_meta('GSE148535')

skin.bk.meta

tdb.sk <- skin.bk |>
  pivot_longer(-1, names_to = 'geo_accession') |>
  left_join(skin.bk.meta) |>
  filter(!is.na(title)) |>
  select(-2) |>
  mutate(group = str_extract(title, '\\d+_hour')) |>
  calc_tpm(group = title, abundance = value) |>
  tidybulk(.sample = title,
           .transcript = Symbol,
           .abundance = value)

tdb.sk |>
  filter(str_detect(title, 'control|162273')) |>
  select(-total.count) |>
  write_csv('mission/FPP/uvb/GSE148535.human.skin.bulk.csv')

tdb.sk <- read_csv('mission/FPP/uvb/GSE148535.human.skin.bulk.csv') |>
  tidybulk(.sample = title,
           .transcript = Symbol,
           .abundance = value)

tdb.sk <- tdb.sk |>
  quick_process_bulk(group = group, skip_scale = T)

pca.sk <- tdb.sk |>
  plot_qc_bulk(value_scaled, group)

pca.sk |>
  filter(!str_detect(title, 'control')) |>
  ggplot(aes(PC1, PC2, color = group, label = title)) +
  geom_point() +
  geom_text_repel() +
  labs(title = 'Human skin UVB bulk sample')

res.sk24 <- tdb.sk |>
  test_differential_abundance(~ 0 + group,
                              contrasts = 'group24_hour-group0_hour',
                              omit_contrast_in_colnames = T)

res.sk6 <- tdb.sk |>
  test_differential_abundance(~ 0 + group,
                              contrasts = 'group6_hour-group0_hour',
                              omit_contrast_in_colnames = T)

res.sk24 |>
  pivot_transcript() |>
  filter(Symbol %in% kegg_mva, FDR < .05)

shorten.sample <- tdb.sk |>
  distinct(title) |>
  mutate(sample = str_c('sample',seq_along(title)))

## export MVA & cytokine ---------
tdb.sk |>
  filter(Symbol %in% c('CCL20','IL6')) |>
  mutate(log.tpm = log1p(tpm)) |>
  select(title, group, Symbol, log.tpm) |>
  pivot_wider(values_from = log.tpm, names_from = Symbol) |>
  write_csv('mission/FPP/uvb/uvb.cytokine.logtpm.csv')

tdb.sk |>
  filter(Symbol %in% kegg_mva) |>
  mutate(zscore = scale(value_scaled)[,1], .by = Symbol) |>
  select(title, group, Symbol, zscore) |>
  pivot_wider(values_from = zscore, names_from = Symbol) |>
  write_csv('mission/FPP/uvb/uvb.MVA.heatmap.csv')

## UPR chaperone ---------
upr.chap <- read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

tdb.sk |>
  filter(Symbol %in% upr.chap$SYMBOL,
         str_detect(title, 'control|162273')) |>
  mutate(group = fct_relevel(group,'0_hour','6_hour')) |>
  as_tibble() |>
  group_by(group) |>
  dplyr::rename(z_score = tpm, gene = Symbol) |>
  heatmap(gene, title, z_score, scale = 'row',
          column_title = 'RNA expression of human skin after UVB treatment',
          palette_value = c('blue','white','red'),
          palette_grouping = list(c('blue','orange','red')),
          show_column_names = F)

up.6h <- res.sk6 |>
  pivot_transcript() |>
  write_csv('mission/FPP/uvb/human.skin.bulk.6h.up.deg.csv')

res.sk6 <- read_csv('mission/FPP/uvb/human.skin.bulk.6h.up.deg.csv')

res.sk6 |>
  filter(Symbol %in% sup.6h$Symbol,
         str_detect(title, 'control'),
         group != '24_hour') |>
  mutate(z.score = scale(value_scaled)[,1], .by = Symbol) |>
  select(title, Symbol, z.score) |>
  pivot_wider(names_from = title, values_from = z.score) |>
  write_csv('uvb.6h.UPR.chaperone.zscore.csv')

sup.6h <- up.6h |>
  filter(Symbol %in% upr.chap$SYMBOL, FDR < .05, logFC > 0)

library(clusterProfiler)
up6h.gsego <- res.sk6 |>
  pull(logFC, name = Symbol) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL',OrgDb = 'org.Hs.eg.db',keyType = 'SYMBOL',
        pvalueCutoff = 1)

uvb_skin_upr_nes <- up6h.gsego@result |>
  as_tibble() |>
  filter(str_detect(Description, 'unfold')) |>
  select(ID, Description, NES, qvalue, leading_edge)

## MVA ---------
tdb.sk |>
  filter(Symbol %in% kegg_mva) |>
  mutate(group = fct_relevel(group,'0_hour','6_hour')) |>
  as_tibble() |>
  group_by(group) |>
  dplyr::rename(z_score = tpm, gene = Symbol) |>
  heatmap(gene, title, z_score, scale = 'row',
          column_title = 'RNA expression of human skin after UVB treatment',
          palette_value = c('blue','white','red'),
          palette_grouping = list(c('blue','orange','red')))

## IFNL -----------
tdb.sk |>
  filter(str_detect(Symbol, '^IFNL')) |>
  mutate(group = fct_relevel(group,'0_hour','6_hour')) |>
  as_tibble() |>
  group_by(group) |>
  rename(z_score = tpm, gene = Symbol) |>
  left_join(shorten.sample) |>
  heatmap(gene, sample, z_score, scale = 'row',
          column_title = 'RNA expression of human skin after UVB treatment',
          palette_value = c('blue','white','red'),
          palette_grouping = list(c('blue','orange','red')))

res.sk24 |>
  pivot_transcript() |>
  filter(str_detect(Symbol, '^IFNL'))

res.sk6 |>
  as_tibble() |>
  filter(str_starts(Symbol, 'IFN'))

## Type I IFN (All IFN but IFNG) -------
tdb.sk |>
  filter(str_detect(Symbol, '^IFN')) |>
  mutate(count_sum = sum(value), .by = Symbol) |>
  filter(count_sum > 0) |>
  mutate(group = fct_relevel(group,'0_hour','6_hour')) |>
  as_tibble() |>
  group_by(group) |>
  dplyr::rename(z_score = tpm, gene = Symbol) |>
  heatmap(gene, title, z_score, scale = 'row',
          column_title = 'RNA expression of human skin after UVB treatment',
          palette_value = c('blue','white','red'),
          palette_grouping = list(c('blue','orange','red')),
          cluster_rows = F,
          show_column_names = F)

## cytokine --------------
tdb.cyt <- tdb.sk |>
  filter(Symbol %in% c('CCL20','IL6'))

sk24.cyt.cusp <- res.sk24 |>
  pivot_transcript() |>
  filter(Symbol %in% c('CCL20','IL6')) |>
  select(Symbol, FDR) |>
  right_join(tdb.cyt) |>
  customize_pvalue(y = log1p(tpm), facets = Symbol,
                   p.value = FDR, group.1 = '0_hour', group.2 = '24_hour')

sk.cyt.cusp <- res.sk6 |>
  pivot_transcript() |>
  filter(Symbol %in% c('CCL20','IL6')) |>
  select(Symbol, FDR) |>
  right_join(tdb.cyt) |>
  customize_pvalue(y = log1p(tpm) + .4, facets = Symbol,
                   p.value = FDR, group.1 = '0_hour', group.2 = '6_hour') |>
  bind_rows(sk24.cyt.cusp)

tpm.sk |>
  filter(Symbol %in% c('CCL20','IL6')) |>
  mutate(group = fct_relevel(group,'0_hour','6_hour')) |>
  ggplot(aes(group, log1p(tpm), color = group)) +
  stat_summary(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~Symbol, scales = 'free_y') +
  theme_pubr() +
  scale_color_manual(values = c('blue','orange','red')) +
  add_manual_pvalue(sk.cyt.cusp) +
  scale_y_continuous(expand = expansion(mult = c(.1,.1))) +
  labs(y = 'log(tpm+1)',
       title = 'RNA expression of human skin after UVB treatment',
       subtitle = 'GSE148535')

last_plot() + theme_jpub

publish_pdf('human.skin.uvb.cytokine.pdf', width = 70)

## lisa2 TF -----
res.sk6 |>
  filter(logFC > 1) |>
  select(Symbol) |>
  write_csv('skin_uvb_6h_up.txt', col_names = FALSE)

res.sk6 |>
  filter(logFC < -1) |>
  select(Symbol) |>
  write_csv('skin_uvb_6h_down.txt', col_names = FALSE)

lisa6up <- read_tsv('basic_motif_skin_uvb_6h_up.txt.lisa.tsv')

lisa6down <- read_tsv('basic_motif_skin_uvb_6h_down.txt.lisa.tsv')

lisa6up |> inner_join(lisa6down,
                     join_by(sample_id, factor)) |>
  mutate(highlight = ifelse(factor %in% c('ATF4','ATF6','XBP1','DDIT3'), factor, NA)) |>
  na.omit() |>
  ggplot(aes(-log10(summary_p_value.x), -log10(summary_p_value.y),
             label = factor, color = highlight)) +
  geom_point() +
  theme_bw()


# primary kera GSE85443 ------
prk43 <- download_ncbi_counts('GSE85443')

prk43.meta <- download_ncbi_meta('GSE85443')

prk43.meta <- prk43.meta |>
  filter(!str_detect(title, '4h')) |>
  mutate(group = str_extract(title, '\\dd...'),
         group = ifelse(is.na(group), 'ctrl', group) |> make.names() |>
           str_replace('X','UV.'),
         sample = str_c(str_extract(title, 'N.'), '_',group),
         title = NULL)

prk43.tdb <- prk43 |>
  select(-2) |>
  pivot_longer(-1, names_to = 'geo_accession') |>
  left_join(prk43.meta) |>
  filter(!is.na(group)) |>
  select(-2) |>
  tidybulk(.sample = sample,
           .transcript = Symbol,
           .abundance = value)

prk43.tdb <- prk43.tdb |>
  preproc_bulk(group)

prk43.tdb |>
  plot_qc_bulk(value_scaled, group)

res.prk43.3d30 <- prk43.tdb |>
  test_differential_abundance(~ 0 + group,
                              contrasts = 'groupUV.3d.30-groupctrl',
                              omit_contrast_in_colnames = T)

res.prk43.3d30 |>
  pivot_transcript() |>
  filter(Symbol %in% key_cytokine, FDR < .05)

prk43.tdb |>
  filter(Symbol %in% key_cytokine) |>
  ggplot(aes(group, value_scaled, color = group)) +
  geom_boxplot() +
  facet_wrap(~Symbol, scales = 'free_y') +
  theme_pubr(x.text.angle = 90, legend = 'right') +
  labs(y = 'Normalized expression',
       title = 'RNA expression of human KC after UVB treatment',
       subtitle = 'GSE85443')

prk43.tdb |>
  filter(Symbol %in% kegg_mva) |>
  ggplot(aes(group, value_scaled, color = group)) +
  geom_boxplot() +
  facet_wrap(~Symbol, scales = 'free_y') +
  theme_pubr(x.text.angle = 90, legend = 'right') +
  labs(y = 'Normalized expression',
       title = 'RNA expression of human KC after UVB treatment',
       subtitle = 'GSE85443')

# mouse skin mck22 -------
mck22 <- read_delim('mission/FPP/uvb/GSE160477_matrix.csv.gz')

mck22

mck22.tdb <- mck22 |>
  pivot_longer(-1) |>
  mutate(group = str_extract(name, 'DETC|KC|LC'),
         time = str_extract(name, '\\d+hr')) |>
  ensembl_to_symbol(.ensembl = ...1) |>
  mutate(transcript = ifelse(is.na(transcript), ...1, transcript)) |>
  select(-ref_genome, -...1) |>
  tidybulk(.sample = name, .transcript = transcript,
           .abundance = value)

mck22.tdb %<>% preproc_bulk(group = c(group,time))

mck22.tdb |> plot_qc_bulk(scaled_abundance = value_scaled, group = group)

res.mck22.tdb <- mck22.tdb |>
  filter(group == 'KC') |>
  test_differential_abundance(~ 0 + time,
                              contrasts = 'time24hr-time0hr',
                              omit_contrast_in_colnames = T)

res.mck22.tdb |>
  pivot_transcript() |>
  filter(transcript %in% str_to_title(c(key_cytokine, kegg_mva)), FDR < .05)

mck22.tdb |>
  filter(transcript %in% str_to_title(kegg_mva), group == 'KC') |>
  mutate(time = fct_relevel(time, '0hr','4hr')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  facet_wrap(~transcript, scales = 'free_y') +
  theme_pubr(x.text.angle = 90, legend = 'right') +
  scale_color_brewer(palette = 'Reds') +
  labs(y = 'Normalized RNA expression',
       title = 'MVA pathway of mouse skin KC after after 300 mJ/cm2 UVB treatment',
       subtitle = 'GSE160477')

mck22.tdb |>
  filter(transcript %in% str_to_title(key_cytokine), group == 'KC') |>
  mutate(time = fct_relevel(time, '0hr','4hr')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  facet_wrap(~transcript, scales = 'free_y') +
  theme_pubr(x.text.angle = 90, legend = 'right') +
  scale_color_brewer(palette = 'Reds') +
  labs(y = 'Normalized RNA expression',
       title = 'Key cytokine of mouse skin KC after 300 mJ/cm2 UVB treatment',
       subtitle = 'GSE160477')

# siege23 mouse 24h -------
sie23 <- readxl::read_excel('mission/FPP/uvb/GSE226792_processed_file_mouse.xlsx')

sie23.key <- sie23 |> select(gene_id, matches('wt.+FPKM')) |>
  pivot_longer(-1) |>
  mutate(group = ifelse(str_detect(name, 'no'), 'ctrl', 'UV'))

sie23.key |>
  filter(gene_id %in% str_to_title(c(key_cytokine))) |>
  ggplot(aes(group, value, color = group)) +
  geom_boxplot() +
  facet_wrap(~gene_id, scales = 'free_y') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  scale_color_hue(direction = -1) +
  labs(y = 'Normalized RNA expression',
       title = 'Key cytokine of mouse epidermis after 24h 100mJ/cm2 UVB treatment',
       subtitle = 'GSE226792')

sie23.key |>
  filter(gene_id %in% str_to_title(c(kegg_mva))) |>
  ggplot(aes(group, value, color = group)) +
  geom_boxplot() +
  facet_wrap(~gene_id, scales = 'free_y') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  scale_color_hue(direction = -1) +
  labs(y = 'Normalized RNA expression',
       title = 'MVA pathway of mouse skin epidermis after 24h 100mJ/cm2 UVB treatment',
       subtitle = 'GSE226792')

sie23.key |>
  filter(gene_id %in% str_subset(sie23.key$gene_id, '^Ifnl')) |>
  ggplot(aes(group, value, color = group)) +
  geom_boxplot() +
  facet_wrap(~gene_id, scales = 'free_y') +
  theme_pubr(x.text.angle = 45, legend = 'right') +
  scale_color_hue(direction = -1) +
  labs(y = 'Normalized RNA expression',
       title = 'IFNL/IFNLR of mouse skin epidermis after 24h 100mJ/cm2 UVB treatment',
       subtitle = 'GSE226792')

sie23.nct <- sie23 |> select(gene_id, matches('wt.+count')) |>
  pivot_longer(-1) |>
  mutate(group = ifelse(str_detect(name, 'noUV'), 'ctrl', 'UV'))

sie23.tdb <- sie23.nct |>
  tidybulk(.sample = name, .transcript = gene_id, .abundance = value)

sie23.tdb %<>% preproc_bulk(group)

sie23.tdb |> plot_qc_bulk(scaled_abundance = value_scaled, group)

sie23.deg <- sie23.tdb |>
  test_differential_abundance(~ 0 + group,
                              contrasts = 'groupUV-groupctrl',
                              omit_contrast_in_colnames = T)

sie23.deg |>
  pivot_transcript() |>
  filter(gene_id %in% str_to_title(c(kegg_mva, key_cytokine)), FDR < .1)

# HaCaT GSE201850 --------
hacat.15 <- read_delim('mission/FPP/uvb/GSE201850.top.table.tsv')

upr.chap <- read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

hacat.15 |>
  filter(Symbol %in% upr.chap$SYMBOL, log2FoldChange < .8) |>
  mutate(log2FoldChange = -log2FoldChange) |>
  ggplot(aes(Symbol, y = '', color = log2FoldChange, size = -log10(padj))) +
  geom_point() +
  theme_pubr() +
  scale_color_distiller(palette = 'RdYlBu') +
  scale_size(range = c(0,3)) +
  theme_jpub +
  labs(title = 'UPR chaperone change in 15 mJ/cm2 UVB treated HaCaT',
       subtitle = 'GSE201850 (n=6)', y = '', x = 'Gene') +
  rotate_x_text(45)

publish_pdf('mission/FPP/figures/hacat.15dose.chaperone.pdf', width = 70)

hacat.15 |>
  filter(Symbol %in% upr.chap$SYMBOL, log2FoldChange < .8) |>
  write_csv('hacat.15dose.chaperone.log2fc.csv')
